No evidence of sustained nonzoonotic Plasmodium knowlesi transmission in Malaysia from modelling malaria case data

Reported incidence of the zoonotic malaria Plasmodium knowlesi has markedly increased across Southeast Asia and threatens malaria elimination. Nonzoonotic transmission of P. knowlesi has been experimentally demonstrated, but it remains unknown whether nonzoonotic transmission is contributing to increases in P. knowlesi cases. Here, we adapt model-based inference methods to estimate RC, individual case reproductive numbers, for P. knowlesi, P. falciparum and P. vivax human cases in Malaysia from 2012–2020 (n = 32,635). Best fitting models for P. knowlesi showed subcritical transmission (RC < 1) consistent with a large reservoir of unobserved infection sources, indicating P. knowlesi remains a primarily zoonotic infection. In contrast, sustained transmission (RC > 1) was estimated historically for P. falciparum and P. vivax, with declines in RC estimates observed over time consistent with local elimination. Together, this suggests sustained nonzoonotic P. knowlesi transmission is highly unlikely and that new approaches are urgently needed to control spillover risks.

Plasmodium knowlesi is a zoonotic malaria parasite carried by simian reservoirs and transmitted to people through bites of infected Anopheles mosquitoes. Since the first report of hundreds of human P. knowlesi cases in Malaysian Borneo in 2004, reported P. knowlesi incidence has markedly increased across Southeast Asia 1,2 . Although increases in cases are partially due to improved detection with the application of molecular diagnostics, data indicate genuine increases in both the proportion and incidence of P. knowlesi 3,4 . In Malaysia, no indigenous cases of nonzoonotic malaria have been reported since 2018 and P. knowlesi is now the main cause of malaria in humans and can cause severe and fatal malaria when untreated 5,6 . These increases in reported incidence are strongly associated with deforestation, suggesting that landscape and land use change may be increasing contact between people, mosquitoes and macaque reservoirs of P. knowlesi 7,8 . P. knowlesi risk in humans is highly spatially heterogenous, with clusters of cases reported in households and villages [9][10][11] . This may reflect similar spatial variation in the force of infection from wildlife reservoir populations or increased adaptation to transmission between human hosts.
Many emerging pathogens have both zoonotic and nonzoonotic transmission pathways in which an initial spillover event from an animal reservoir may lead to stuttering chains of human-to-human transmission or larger outbreaks 12 . For P. knowlesi, genetic evidence suggests that transmission to humans remains primarily zoonotic, occurring through bites from vectors infected by wild non-human primate populations 13 . Evidence suggests P. knowlesi infections are typically chronic with lowlevel parasitemias in long-tailed macaques (Macaca fascicularis), the main reservoir 14 . The spatial distribution of human cases mirrors P. knowlesi infection rates in macaques, with the highest infection rates in humans and macaques reported in Malaysian Borneo 15 . P. knowlesi has a rapid 24 h replication rate with observed prepatent periods in humans ranging from 9 to 12 days 16 . Although limited data are available on the duration and frequency of human infectiousness to mosquitoes, the feasibility of human-mosquito-human transmission of P. knowlesi has been demonstrated once experimentally in a laboratory study 16 but has not been formally described in natural settings. A recent systematic literature review concluded nonzoonotic P. knowlesi transmission was biologically plausible, but there was no empirical evidence that humanmosquito-human P. knowlesi transmission occurs naturally 17 . Notably, no previous studies had systematically assessed whether spatiotemporal patterns of P. knowlesi cases were consistent with nonzoonotic transmission.
Quantifying the relative contributions of zoonotic and nonzoonotic transmission is essential for designing effective surveillance and control. Transmission pathways are represented by two epidemiological parameters: the spillover rate, the rate at which a pathogen is transmitted from animals to humans, and the reproductive number (R), the number of secondary human cases resulting from an infectious human individual over the course of their infection 18 . Directly attributing human infections to zoonotic sources is challenging as the distributions and prevalence of wildlife reservoirs are largely unknown, and infection occurs indirectly via mosquito vectors. Alternatively, model-based inference methods can be used to infer likely transmission routes from routinely collected surveillance data on reported cases [18][19][20][21] . Diffusion network approaches have been used to quantify individual malaria case reproduction numbers (R C ) by evaluating the probability of two cases being connected using the time of symptom onset, spatial locations of cases, and estimated serial interval (SI), i.e. the time from symptom onset in one case to symptom onset in the secondary case 19,[22][23][24] . These methods also allow assessment of the likelihood of cases occurring from outside (unobserved) sources not represented in human surveillance data, such as introductions through zoonotic spillover 18,21 .
Here, we use a diffusion network approach within a Bayesian framework to estimate R C for P. knowlesi and nonzoonotic malaria species across Malaysia from 2012 to 2020. Due to an active malaria elimination programme, Malaysia has an exceptionally strong surveillance system, and detailed information is available for all cases. We use this surveillance dataset and parameters identified from the literature to estimate the SI for nonzoonotic P. knowlesi transmission. We then quantify estimates and uncertainty around R C for the zoonotic malaria P. knowlesi and nonzoonotic malaria P. falciparum and P. vivax. Models of P. knowlesi transmission were consistent with a large reservoir of unobserved infection sources, suggesting human cases are primarily driven by spillover and with limited subcritical onward transmission (R C < 1). In contrast, R C estimates for P. falciparum and P. vivax identified historical sustained human-mosquito-human transmission (i.e. R C > 1), with subsequent decreases in R C observed over time as the country approached elimination of transmission of these two nonzoonotic malaria parasites.

Results
Between 16 December 2011 and 3 January 2021, 32,635 malaria cases were reported to the Malaysian national malaria surveillance system. Of these, the majority (n = 26,093, 79.95%) were reported from the East Malaysian states of Sabah and Sarawak (Fig. 1). All suspected P. knowlesi cases were confirmed by molecular methods, with confirmed P. knowlesi accounting for over 70% (n = 23,143) of malaria cases nationally and no indigenous nonzoonotic malaria cases reported after 2018 (Fig. 2). Based on travel history to malaria-endemic areas within the past month and previous malaria diagnoses, cases are classified by the Malaysian Ministry of Health as indigenous, introduced, imported or reoccurrences (for P. vivax). Most malaria cases were classified as indigenous (n = 27,125/32,635, 83.12%), including over 99% of P. knowlesi cases. Over 80% (n = 18,691/23,143) of P. knowlesi cases were reported in men, with a median age of 38 years (interquartile range [IQR]: 26-49 years). Similarly, most nonzoonotic malaria cases were also reported in adult men; the median age was 31 (IQR: 21-43 years), and 84% (7945/9489) were male.
Malaysia has high health system coverage with free treatment for malaria 25 . For available surveillance data, the median time between symptom onset and treatment was 5 days, with over 97% of individuals seeking treatment within 2 weeks (Fig. 3a). Out of 16,765 P. knowlesi cases screened for gametocytes by microscopy, 5.74% (95% confidence interval [CI]: 5.39-6.10%) had microscopically detectable gametocytes, consistent with other reports on the P. knowlesi gametocyte carriage in clinical cases 26 . Infectiousness to mosquitoes is primarily dependent on the density of malaria gametocytes, with infectious individuals predominantly having gametocyte densities high enough to be detected microscopically 27,28 . Proportions of infectious P. knowlesi cases increased with time since symptom onset; however, due to the relatively rapid reporting and treatment times, the mean duration of human infectiousness was only 1.645 days (range: 0-27 days) (Fig. 3b, Supplementary Information).
We assessed the likelihood of cases being part of the same transmission chain based on species-specific SI parameters, dates and geographical locations of reported cases, and case classification as indigenous or non-indigenous (imported, introduced or    reoccurrences). This modelling framework additionally included a threshold parameter (ε) that prevents the likelihood from becoming zero due to a complete break in the transmission chain and, therefore, epidemiologically represents the probability of infection from unobserved external sources needed to sustain transmission. We assessed ε ranges using mean prior values ranging from <0.001, a prior distribution centring around the assumption that nearly all cases are infected from human sources included in the surveillance data to 1, a prior distribution assuming that all human cases are infected from external zoonotic or imported sources (Supplementary Table 4) 19 . P. knowlesi models with ε values approaching 1 resulted in better fits, implying that all or nearly all human cases are infected from external sources, most likely to be zoonotic spillover from non-human primates in this context (Supplementary Tables 4  and 5). For each individual, their estimated number of secondary infections was calculated to give a time-varying point estimate of R C . An individual R C value < 1 indicates a low likelihood that the case is the source infection for any other case in the dataset. Mean R C values for P. knowlesi in East Malaysia from 2012 to 2020 were 0.0739 (95% CI: 0.0006-0.278) and 0.0147 (95% CI: 0.000003-0.0637) for West Malaysia. There were no estimates of R C greater than 1 in the bestfitting models, and no scenarios assessed that had mean R C values greater than 1, indicating no sustained transmission (Supplementary Information). As we estimated time-varying R C values for each individual, the estimates of all R C values as less than 1 also suggest there is no major overdispersion, excluding large super-spreading events.
Although there were clear increasing trends in the numbers of reported P. knowlesi cases through the study period, this was not reflected in increases in R C estimates, and there was no evidence of deviation from zoonotic-mosquito-human transmission patterns.
In contrast, models for nonzoonotic malaria identified evidence for clear indigenous human-mosquito-human transmission chains, with 17% of P. falciparum cases and 24% of P. vivax cases with estimated R C values greater than 1 from 2012 to 2020. For P. falciparum, the mean R C value was 0.568 (95% CI: 0-1.97) for East Malaysia and 0.386 (95% CI: 0-2.04) for West Malaysia. The mean R C for P. vivax was 0.524 (95% CI: 0-2.436) and 0.230 (95%CI: 0, 1.820) for East and West Malaysia, respectively (Fig. 4). The probability of P. falciparum or P. vivax cases leading to onward transmission varied over space and time (Fig. 5). There were clear decreasing temporal trends in R C estimates for both nonzoonotic malaria species, with no R C values above 1 detected after 2018 when the last indigenous nonzoonotic malaria case was reported by the Malaysian Ministry of Health. Models clearly identified an outbreak of P. vivax reported in West Malaysia between 2016 and 2017, demonstrating the utility of these methods in pre-elimination settings 29 .

Discussion
Quantifying rates of nonzoonotic transmission is critical to understanding the dynamics of emerging zoonotic diseases and identifying appropriate interventions. While many vector-borne and zoonotic diseases have emerged to become widespread in human populations (e.g. P. knowlesi, Lyme disease), the increase in human case numbers is not necessarily correlated with increased nonzoonotic transmission 30,31 . Here, we demonstrate the use of quantitative modelbased approaches to identify key parameters and estimate individual reproductive numbers of the zoonotic and vector-borne disease P. knowlesi. Model results strongly suggest that P. knowlesi remains a primarily zoonotic disease driven by spillover, finding no evidence of sustained nonzoonotic P. knowlesi transmission. These results sharply contrast with R C estimates of the nonzoonotic malaria P. falciparum and P. vivax, where clear chains of human-mosquito-human transmission were identified. This highlights the utility of this method to assess the probability of nonzoonotic transmission using high-quality routine surveillance data and provides a critical tool to monitor changing transmission patterns.
A key determinant of the probability of nonzoonotic P. knowlesi transmission is whether humans are infectious. Using the largest dataset of gametocyte carriage rates in P. knowlesi clinical cases reported, consistent proportions of cases had microscopically observed gametocytes (5.74%), supporting the possibility of nonzoonotic P. knowlesi transmission. However, P. knowlesi gametocyte carriage rates in clinical cases upon presentation to health facilities are notably lower than for other nonzoonotic malaria species. For example, for P. vivax infections, which are biologically similar to P. knowlesi, up to 100% of clinical cases have microscopically observed Malaysia. Data are presented as the median, first and third quartiles and range, with points more than 1.5 times the interquartile range plotted separately as outliers.
Article https://doi.org/10.1038/s41467-023-38476-8 gametocytes at presentation 32 . More extensive meta-analyses of gametocyte carriage in P. falciparum cases have also identified a wider range (12.1-32.2%) of microscopically observed gametocytes at rates typically higher than P. knowlesi cases reported within this dataset 33 . While this does suggest human P. knowlesi cases have the potential to be infectious, this parasite may be less adapted to human hosts in terms of gametocytogenesis. Critical questions remain about how infective P. knowlesi-infected humans are to mosquitoes, particularly considering the strong circadian rhythms of infectiousness to mosquitoes observed in P. knowlesi-infected non-human primates 34 . Our model fits strongly suggests that P. knowlesi remains primarily driven by zoonotic spillover, finding no evidence of sustained nonzoonotic P. knowlesi transmission. For all of the best fitting models, there were no estimates of R C greater than 1, and in all scenarios assessed, there were no estimates suggesting R C is routinely over 1. There was additionally no clear evidence of changes to P. knowlesi R C over time that would suggest P. knowlesi may be adapting to transmission between humans. This contrasts with R C estimates for both P. falciparum and P. vivax, where R C estimates frequently exceed 1 and show clear decreasing temporal trends during an extensive malaria elimination campaign. This highlights the utility of these methods in monitoring changing transmission patterns of both emerging and endemic vector-borne diseases 19,20 .
While these results can shed insight into broad transmission patterns, this analysis is not without significant limitations. Modelbased estimates alone cannot be used to conclusively prove that two cases are linked. This modelling approach estimated R C under different assumptions but did not reconstruct transmission networks linking cases. Further genetic, epidemiological and ecological data would be needed to prove transmission occurred between two individuals and rule out multiple spillover events occurring from zoonotic reservoirs. Additionally, while P. knowlesi results are consistent with a large reservoir of unobserved infections, this modelling framework does not explicitly attribute these to zoonotic sources. Other potential sources include asymptomatic or unreported infections. However, given the high coverage of the Malaysian surveillance system, infrequency and low parasite density of asymptomatic human P. knowlesi infections 35,36 and widespread distribution of infected simian reservoirs 15 , unobserved infections would appear to be mostly from zoonotic sources. As higher rates of low-density, asymptomatic P. knowlesi infections have been reported in people residing near clinical cases, further studies are needed to identify whether peri-domestic or localised nonzoonotic P.  The probability of R C exceeding 1 is shown for a P. falciparum and b P. vivax; P. knowlesi probabilities were not calculated as there were no estimates of R C greater than 1 in the best-fitting models, but village-level R C estimates can be found in Supplementary Fig. 3.
knowlesi transmission is occurring but not detected by clinical surveillance systems 35 . Higher estimates of R C values could be used to guide targeted investigations at locations with the highest likelihood of nonzoonotic transmission. Future studies could additionally explore geographic variability in human-mosquito contact rates, human immunity, treatment-seeking behaviours and other factors impacting transmission probabilities and SI distributions. Despite these limitations, results indicate there is no evidence of sustained nonzoonotic P. knowlesi transmission and demonstrate decreases in the transmission of P. falciparum and P. vivax in Malaysia throughout the past decade. This has important implications for the design of malaria surveillance and control programmes in areas with zoonotic malaria. Current measures, such as insecticide-treated net distribution and active case detection, have proven highly effective at controlling nonzoonotic malaria. While these measures may also be reducing nonzoonotic P. knowlesi transmission, new control strategies are urgently needed to mitigate the risks of spillover from wildlife hosts. With the inclusion of zoonotic malaria in global malaria elimination certification criteria, P. knowlesi poses a major threat to malaria elimination across Southeast Asia 37 . These methods provide a template for monitoring potential changes in transmission patterns of P. knowlesi and other zoonotic and vector-borne diseases.

Surveillance data
National malaria surveillance data collected by the Malaysian Ministry of Health were assembled from 16 December 2011 to 3 January 2021. Data on malaria cases included demographic data (age, gender and occupation), case notification date to the Ministry of Health, symptom onset date, diagnosis date, method of detection, importation status, parasite species diagnosed, presence of gametocytes and health facility of detection. Symptom onset dates were excluded as unreliable if they occurred after diagnosis or notification dates or more than 30 days prior to diagnosis. Due to frequent misdiagnosis of P. knowlesi by microscopy, diagnosis of all presumptive P. knowlesi cases and imported malaria, cases were confirmed using species-specific molecular methods 25 . Additionally, surveillance data included geographical data on the case residence, probable location of the infection (state, district, village name and/or GPS coordinates) and travel history. All case locations were manually confirmed and geolocated to the nearest village centroid (Supplementary Information). Cases without accurate address data were excluded from the final analysis (n = 2117).

SI estimation
The SI, the time between symptom onset of primary and secondary malaria cases, requires estimating the duration of a series of sequential processes which need to occur in a human-mosquito-human transmission cycle 38 . While these intervals can be inferred from contract tracing data for directly transmitted diseases, SIs can only be estimated indirectly for vector-borne diseases 39,40 . SI durations have previously been estimated for P. falciparum and P. vivax 19,40 . However, there is a lack of empirical evidence on human-mosquito-human P. knowlesi transmission as this has only been observed once under experimental conditions 16 . The rapid replication cycle of P. knowlesi and weak evidence of adaptation to humans suggests the SI may differ from that of other nonzoonotic malaria parasites 34 .
To estimate the SI for human-mosquito-human P. knowlesi transmission, we adapted a quantitative model-based approach 40 . This models the SI as the sum of random variables representing sequential steps in the transmission cycle, including the prepatent period, human to mosquito transmission period, extrinsic incubation period, mosquito-to-human transmission period and infection-to-reporting periods. These were parameterised based on a combination of data from secondary literature and the Malaysian malaria surveillance dataset (Supplementary Information). Due to the typically low parasite densities and infrequency of asymptomatic human P. knowlesi infections in the region, we assumed asymptomatic infections did not contribute to transmission 35,41 . We also assumed individuals were not infectious after reporting due to Malaysia's mandatory treatment and hospitalisation policy. Finally, no data were available on the duration or timing of infectiousness of P. knowlesi in humans, and the relatively short duration between symptoms onset and treatment within this population precluded estimation of natural recovery rates. While P. knowlesi gametocytes are highly synchronous and have strong circadian rhythms in macaques 34,42 , available evidence suggests this is not the case in human infections 43 .

Transmission likelihood and R C estimation
To estimate the joint likelihood of transmission between cases and a case having an unobserved source of infection, we used an adapted version of the NetRate algorithm to calculate R C for each case 19,20,22,44 . This uses a Bayesian framework to estimate the likelihood of two cases being within the transmission chain based on the geographical locations of cases, indigenous case status, time of symptom onset and SI parameters (Supplementary Information). Case reproduction numbers (R C ) are then estimated by, for each case, summing the normalised likelihood of that case being the infector of each of the other cases. Non-indigenous cases and cases with a symptom onset date before the candidate case are not included as potential infectees. This approach additionally includes a provision for missing or unobserved sources of infections that can infect any observed individual with probability ε. These epsilon edges (ε) link cases to an external source when the likelihood of transmission between two cases is sufficiently low. Higher values of ε indicate a larger proportion of cases are assumed to be infected by external sources outside the surveillance data; in the case of P. knowlesi, this is mostly to be zoonotic reservoirs.
We fit separate models for P. knowlesi, P. falciparum and P. vivax for both East Malaysia (Borneo) and West Malaysia. From surveillance reports, most cases reported only local travel within the same village or subdistrict. To account for a high likelihood of local transmission, we used a fixed value for δ of 0.1, the parameter determining transmission likelihood based on Euclidean distance between cases (assessed in kilometres); this corresponds to most transmission events occurring within 10 km based on reported travel history (Supplementary Information). We fit shifted Rayleigh distributions to describe a prior distribution of possible SIs for P. knowlesi and nonzoonotic malaria. This distribution is defined by α, the shape of the SI distribution, and γ, the shifting parameter accounting for the minimum time between a bite from an infected mosquito and infectiousness in a human. For P. falciparum and P. vivax, we used previously published parameters to inform the SI distribution, with γ fixed at 15 days and normally distributed priors on α with a mean of 0.003 and a standard deviation of 0.1; this corresponds to a mean SI of 36 days (95% range: 21-60 days) 20 . For P. knowlesi, the γ was fixed at 14 days based on the modelled time to infectiousness fit from the Malaysian surveillance data (Supplementary Information). To account for uncertainty in the SI distribution, we used a normally distributed prior for α with a mean of 0.002 and standard deviation of 0.001, capturing the full range of plausible estimates of the P. knowlesi SI.
For each malaria species and location, surveillance data was input as a time-ordered series of malaria notification dates. As there was considerable uncertainty in the proportion of zoonotic transmission of P. knowlesi, we conducted a sensitivity analysis using normally distributed priors for ε with mean values of 0.01, 0.1 and 1. As ε values were bounded between 1e −11 and 1, prior estimates were truncated at 1, forming a half-normal estimate for higher values. Similar sensitivity analysis was conducted for P. falciparum and P. vivax, with mean values of priors for ε ranging from 0.0001 to 0.000001, reflecting the expected high sensitivity of malaria surveillance systems. We assessed model fit based on second-order Akaike's Information Criteria (AIC) and visually inspected the posterior distributions of key parameters (γ, δ, ε) to select the best-fitting model. For all final models, we estimated R C by creating a matrix of potential infectors and infectees with the normalised likelihood of cases being connected. For each individual, the fractional number of secondary infections was calculated to give a time-varying point estimate of R C . Models were previously validated through simulation studies 44 , with code available at https://github. com/IzzyRou/spatial_rcs. Spatial and temporal patterns of R C were explored for all species and locations. For P. falciparum and P. vivax, the probability of onward transmission (R C > 1) was assessed using statistical models with joint temporally and spatially structured random effects with Bayesian inference approximated using the Integrated Nested Laplace Approximation (INLA) 45 . We classified cases into dichotomous categories based on R C estimates. For each species, we fit a binomial model with a logit link to the number of cases with R C estimates > 1 out of the total cases per year detected within 5 km 2 grid cells (Supplementary Information). The best-fitting model was selected using Deviance Information Criteria (DIC) and included a spatial effect and a temporally structured random walk model of order 2 46 . The spatial effect was modelled as a Matern covariance function using an approximation to a stochastic partial differential equations approach. Posterior expectations and probabilities were estimated using 1000 samples. All data analysis was conducted in R statistical software v 4.1, and data were checked and cleaned using Quantum GIS v 3.24.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
As malaria case surveillance datasets contain identifiable health information, data is only available following approval of the Medical Ethics and Research Committee in Malaysia, relevant ethics committees in the UK and with permission from the Malaysian Ministry of Health. Requests to access datasets should be directed to Kimberly.Fornace@lshtm.ac.uk.